## Replication File for:
## ``Concentrated Burdens: How Self-Interest and Partisanship Shape Opinion on Opioid Treatment Policy"
## Justin de Benedictis-Kessner and Michael Hankinson


## ---------- ##
#### Notes: ####
## ---------- ##

# This is the second script file needed to reproduce the results from the paper
# Run the `01-survey-recodes.R`` script first to produce recoded survey dataset that this script needs

## --------------------- ##
#### Preliminary stuff ####
## --------------------- ##

rm(list=ls())
library(stargazer)
library(tidyverse)
library(xtable)

library(Cairo)
# CairoFonts(regular="CMU Serif",
# 					 bold="CMU Serif Extra"
# )

red_mit = '#A31F34'
red_light = '#A9606C'
blue_mit = '#315485'
grey_light= '#C2C0BF'
grey_dark = '#8A8B8C'
black = '#353132'

## ----------------- ##
#### Recoded data: ####
## ----------------- ##

a <- read_csv("norc_data_recoded.csv")

## ----------------------- ##
#### Sample Demographics ####
## ----------------------- ##

## Table A1:
stargazer(as.data.frame(a),summary = T,
					keep = c(
						"college$",
						"married$",
						"income_gt50k$",
						"AGE$",
						"female$",
						"^white$",
						"^black$",
						"^hispanic$",
						"dem",
						"liberal",
						"own",
						"above_income",
						"above_death_rate",
						"personal_bi"),
					covariate.labels = c(
						"Age",
						"\\% Above Median Overdose Rate",
						"\\% Above State Median Income",
						"\\% Female",
						"\\% Democrat",
						"\\% Liberal",
						"% Income >50k",
						"\\% College Degree",
						"\\% Married",
						"\\% White",
						"\\% Black",
						"\\% Hispanic",
						"\\% Homeowner",
						"\\% Know Someone with Addiction"
					),
					float = F,nobs = F,
					out = "descriptives.tex")

fit_pers <- lm(personal_bi ~ AGE + female + black + hispanic + own  + factor(pid) + factor(ideo) + college + married+ income_gt50k + factor(region), data=a)


## Table A2:
stargazer(fit_pers,
					dep.var.labels = "Know Someone w/ Addiction",
					covariate.labels = c("Age",
															 "Female",
															 "Black",
															 "Hispanic",
															 "Homeowner",
															 "Party: Independent",
															 "Party: Republican",
															 "Ideology: Moderate",
															 "Ideology: Conservative",
															 "College Degree",
															 "Married",
															 "Income $>$50k",
															 "Region: Northeast",
															 "Region: South",
															 "Region: West"),
					notes = c("Omitted category for partisanship is 'Democrat',","for ideology is 'Liberal', and for Region is","'Midwest'"),
					notes.align = "l",
					float = F,df=F,omit.stat = c("rsq","ser"),out = "reg_tab_pers.tex"
)

## Figure A1:
personal_hist <- ggplot(subset(a,Q3!=98)) + 
	geom_bar(aes(x=Q3,y=(..count..)/sum(..count..))) + 
	scale_y_continuous("% responding",labels=scales::percent,limits=c(0,0.5)) + 
	scale_x_discrete("Know anyone with opioid addiction?",
									 limits=c(1:5),labels=c("Self","Family\nmember","Close\nfriend","Acquaintance","No one")) + 
	theme_bw() + 
	theme(axis.text = element_text(size=10))

CairoPDF("personal_hist.pdf", width=7, height=5)
print(personal_hist)
dev.off()


#### Balance across treatment conditions: ####
library(cobalt)
a$condition01 <- ifelse(a$distance=="far",1,0)
btab <- bal.tab(covs= a[,c("AGE", "above_death_rate", "above_income", "female", "dem", "rep", "liberal", "conservative", "income_gt50k", "college", "married", "white", "black", "hispanic", "own", "personal_bi")],
								treat="condition01",binary = "raw",stat="mean.diffs",
								data=a)

btab_formatted <- btab$Balance[!grepl("<NA>",rownames(btab$Balance)),c("M.0.Un","M.1.Un")]
rownames(btab_formatted) <- c(
	"Age",
	"% Above Median Overdose Rate",
	"% Above State Median Income",
	"% Female",
	"% Democrat",
	"% Republican",
	"% Liberal",
	"% Conservative",
	"% Income >50k",
	"% College Degree",
	"% Married",
	"% White",
	"% Black",
	"% Hispanic",
	"% Homeowner",
	"% Know Someone w/ Addiction"
)
colnames(btab_formatted) <- c("Mean[Near condition]","Mean[Far condition]")

# t tests:
treat_age <- t.test(x=a$AGE[a$distance=="far"],y=a$AGE[a$distance=="near"])
treat_death <- t.test(x=a$above_death_rate[a$distance=="far"],y=a$above_death_rate[a$distance=="near"])
treat_income <- t.test(x=a$above_income[a$distance=="far"],y=a$above_income[a$distance=="near"])
treat_female <- t.test(x=a$female[a$distance=="far"],y=a$female[a$distance=="near"])
treat_dem <- t.test(x=a$dem[a$distance=="far"],y=a$dem[a$distance=="near"])
treat_rep <- t.test(x=a$rep[a$distance=="far"],y=a$rep[a$distance=="near"])
treat_lib <- t.test(x=a$liberal[a$distance=="far"],y=a$liberal[a$distance=="near"])
treat_cons <- t.test(x=a$conservative[a$distance=="far"],y=a$conservative[a$distance=="near"])
treat_income_gt50k <- t.test(x=a$income_gt50k[a$distance=="far"],y=a$income_gt50k[a$distance=="near"])
treat_college <- t.test(x=a$college[a$distance=="far"],y=a$college[a$distance=="near"])
treat_married <- t.test(x=a$married[a$distance=="far"],y=a$married[a$distance=="near"])
treat_white <- t.test(x=a$white[a$distance=="far"],y=a$white[a$distance=="near"])
treat_black <- t.test(x=a$black[a$distance=="far"],y=a$black[a$distance=="near"])
treat_hisp <- t.test(x=a$hispanic[a$distance=="far"],y=a$hispanic[a$distance=="near"])
treat_own <- t.test(x=a$own[a$distance=="far"],y=a$own[a$distance=="near"])
treat_pers <- t.test(x=a$personal_bi[a$distance=="far"],y=a$personal_bi[a$distance=="near"])

btab_formatted$`p-value of difference` <- c(treat_age$p.value,
																						treat_death$p.value,
																						treat_income$p.value,
																						treat_female$p.value,
																						treat_dem$p.value,
																						treat_rep$p.value,
																						treat_lib$p.value,
																						treat_cons$p.value,
																						treat_income_gt50k$p.value,
																						treat_college$p.value,
																						treat_married$p.value,
																						treat_white$p.value,
																						treat_black$p.value,
																						treat_hisp$p.value,
																						treat_own$p.value,
																						treat_pers$p.value
)

## Output Table A3:
print(xtable(btab_formatted),booktabs=T,file="balance_table.tex",floating = F)


# -------------------------- #
#### Descriptive Analyses ####
# -------------------------- #

# party id: effect of party id
redist_pid <- lm(ability_support_bi ~ factor(pid), data=a);summary(redist_pid)

writeLines(as.character(100*abs(round(redist_pid$coefficients[3],2))), "redist_pid_perc.tex",sep = "")

needs_pid <- lm(needs_support_bi ~ factor(pid), data=a);summary(needs_pid)

writeLines(as.character(100*abs(round(needs_pid$coefficients[3],2))), "needs_pid_perc.tex",sep = "")

nimby_pid <- lm(dist_support_bi ~ factor(pid), data=a);summary(nimby_pid)

writeLines(as.character(100*abs(round(nimby_pid$coefficients["factor(pid)1"],2))), "nimby_pid_perc.tex",sep = "")

# ideology as a robustness check
redist_ideo <- lm(ability_support_bi ~ factor(ideo), data=a);summary(redist_ideo)
needs_ideo <- lm(needs_support_bi ~ factor(ideo), data=a);summary(needs_ideo)

# based on personal exposure:
writeLines(as.character(100*round(mean(a$personal_bi, na.rm=T),2)),"n_pers1_perc.tex",sep = "")

# Output average support for three DVs:
writeLines(as.character(100*round(mean(a$support_bi[a$condition=="ability"],na.rm=T),2)),"average_ability_perc.tex",sep = "")
writeLines(as.character(100*round(mean(a$support_bi[a$condition=="need"],na.rm=T),2)),"average_needs_perc.tex",sep = "")
writeLines(as.character(100*round(mean(a$dist_support_bi,na.rm=T),2)),"average_nimby_perc.tex",sep = "")




#### Figure 1: Descriptive Differences for all three policy DVs
desc_full <- a %>%
	select(needs_support_bi,
				 ability_support_bi,
				 dist_support_bi) %>%
	gather(needs_support_bi,
				 ability_support_bi,
				 dist_support_bi,
				 key = "variable",value = "value") %>%
	group_by(variable) %>%
	summarize(mean = mean(value,na.rm=T),
						sd = sd(value,na.rm=T),
						n = sum(!is.na(value))
	) %>%
	mutate(ci.95.lo = mean + qnorm(0.025)*sd/sqrt(n),
				 ci.95.hi = mean + qnorm(0.975)*sd/sqrt(n),
				 Subgroup = "All\nRespondents"
	)

desc_pid <- a %>%
	filter(!is.na(pid)) %>%
	select(pid,
				 needs_support_bi,
				 ability_support_bi,
				 dist_support_bi) %>%
	gather(needs_support_bi,
				 ability_support_bi,
				 dist_support_bi,
				 key = "variable",value = "value") %>%
	group_by(pid,variable) %>%
	summarize(mean = mean(value,na.rm=T),
						sd = sd(value,na.rm=T),
						n = sum(!is.na(value))
	) %>%
	mutate(ci.95.lo = mean + qnorm(0.025)*sd/sqrt(n),
				 ci.95.hi = mean + qnorm(0.975)*sd/sqrt(n),
				 Subgroup = recode(pid,`-1`="Democrats",`0`="Independents",`1`="Republicans")
	) %>%
	ungroup() %>%
	select(-pid)

desc_ideo <- a %>%
	filter(!is.na(ideo)) %>%
	select(ideo,
				 needs_support_bi,
				 ability_support_bi,
				 dist_support_bi) %>%
	gather(needs_support_bi,
				 ability_support_bi,
				 dist_support_bi,
				 key = "variable",value = "value") %>%
	group_by(ideo,variable) %>%
	summarize(mean = mean(value,na.rm=T),
						sd = sd(value,na.rm=T),
						n = sum(!is.na(value))
	) %>%
	mutate(ci.95.lo = mean + qnorm(0.025)*sd/sqrt(n),
				 ci.95.hi = mean + qnorm(0.975)*sd/sqrt(n),
				 Subgroup = recode(ideo,`-1`="Liberals",`0`="Moderates",`1`="Conservatives")
	) %>%
	ungroup() %>%
	select(-ideo)

desc_persexp <- a %>%
	filter(!is.na(personal_bi)) %>%
	select(personal_bi,
				 needs_support_bi,
				 ability_support_bi,
				 dist_support_bi) %>%
	gather(needs_support_bi,
				 ability_support_bi,
				 dist_support_bi,
				 key = "variable",value = "value") %>%
	group_by(personal_bi,variable) %>%
	summarize(mean = mean(value,na.rm=T),
						sd = sd(value,na.rm=T),
						n = sum(!is.na(value))
	) %>%
	mutate(ci.95.lo = mean + qnorm(0.025)*sd/sqrt(n),
				 ci.95.hi = mean + qnorm(0.975)*sd/sqrt(n),
				 Subgroup = recode(personal_bi,`1`="Know\nSomeone\nw/ Addiction",`0`="Don't Know\nSomeone\nw/ Addiction")
	) %>%
	ungroup() %>%
	select(-personal_bi)



desc_income <- a %>%
	filter(!is.na(above_income)) %>%
	select(above_income,
				 needs_support_bi,
				 ability_support_bi,
				 dist_support_bi) %>%
	gather(needs_support_bi,
				 ability_support_bi,
				 dist_support_bi,
				 key = "variable",value = "value") %>%
	group_by(above_income,variable) %>%
	summarize(mean = mean(value,na.rm=T),
						sd = sd(value,na.rm=T),
						n = sum(!is.na(value))
	) %>%
	mutate(ci.95.lo = mean + qnorm(0.025)*sd/sqrt(n),
				 ci.95.hi = mean + qnorm(0.975)*sd/sqrt(n),
				 Subgroup = recode(above_income,`1`="Above Median\nIncome",`0`="Below Median\nIncome")
	) %>%
	ungroup() %>%
	select(-above_income)

desc_death <- a %>%
	filter(!is.na(above_death_rate)) %>%
	select(above_death_rate,
				 needs_support_bi,
				 ability_support_bi,
				 dist_support_bi) %>%
	gather(needs_support_bi,
				 ability_support_bi,
				 dist_support_bi,
				 key = "variable",value = "value") %>%
	group_by(above_death_rate,variable) %>%
	summarize(mean = mean(value,na.rm=T),
						sd = sd(value,na.rm=T),
						n = sum(!is.na(value))
	) %>%
	mutate(ci.95.lo = mean + qnorm(0.025)*sd/sqrt(n),
				 ci.95.hi = mean + qnorm(0.975)*sd/sqrt(n),
				 Subgroup = recode(above_death_rate,`1`="Above Median\nOverdose Rate",`0`="Below Median\nOverdose Rate")
	) %>%
	ungroup() %>%
	select(-above_death_rate)

## add all summaries together
desc_tab <- bind_rows(desc_full,
											desc_pid,
											desc_income,
											desc_death,
											desc_persexp)

desc_tab$Subgroup <- factor(desc_tab$Subgroup,
														levels = c("All\nRespondents",
																			 "Democrats","Independents","Republicans",
																			 # "Income <$50k","Income >$50k",
																			 # "Liberals","Moderates","Conservatives",
																			 "Know\nSomeone\nw/ Addiction","Don't Know\nSomeone\nw/ Addiction",
																			 
																			 "Below Median\nIncome",
																			 "Above Median\nIncome",
																			 "Below Median\nOverdose Rate",
																			 "Above Median\nOverdose Rate"
														),
														ordered = T)
desc_tab$Policy <- factor(desc_tab$variable,levels = c("ability_support_bi","needs_support_bi","dist_support_bi"),
													labels=c("Treatment Funding: Redistributive","Treatment Funding: Needs-Based","Clinic Construction"),
													ordered=T)
desc_tab$n_text <- paste0("n = ",desc_tab$n)


desc_tab$Subgroup_n <- paste0(desc_tab$Subgroup,"\n(",desc_tab$n_text,")")


# main effects points
x_labels <- desc_tab[which(desc_tab$variable=="ability_support_bi"),c("Subgroup","Subgroup_n")]
plot_desc_points_a <- ggplot(subset(desc_tab,variable=="ability_support_bi"),aes(x=Subgroup,y=mean)) + 
	geom_errorbar(aes(ymin=ci.95.lo, ymax=ci.95.hi),
								position=position_dodge(width=0.5),
								width=0.25, size=0.5) +
	geom_point(aes(shape=Policy), 
						 position=position_dodge(width=0.5),size=1.5) + 
	scale_y_continuous("% Support Policy",limits=c(0.2,0.83),
										 breaks=seq(0,1,0.2),
										 labels=100*seq(0,1,0.2)) + 
	scale_x_discrete("Subgroup",breaks=x_labels$Subgroup,labels=x_labels$Subgroup_n) + 
	facet_wrap("Policy",ncol = 1) + 
	# geom_text(data=desc_tab,aes(x = Subgroup,y=0.20,label=n)) +
	theme_bw() + 
	theme(legend.position = "none",
				strip.text = element_text(size=18),
				axis.text = element_text(size=12))

fname = "descriptive_points_panel_a.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_desc_points_a, width=11, height=5)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=11, height=5)
	print(plot_desc_points_a)
	dev.off()
}

fname = "descriptive_points_panel_a_withpers.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_desc_points_a, width=12.75, height=5)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=12.75, height=5)
	print(plot_desc_points_a)
	dev.off()
}

# panel b
x_labels <- desc_tab[which(desc_tab$variable=="needs_support_bi"),c("Subgroup","Subgroup_n")]
plot_desc_points_b <- ggplot(subset(desc_tab,variable=="needs_support_bi"),aes(x=Subgroup,y=mean)) + 
	geom_errorbar(aes(ymin=ci.95.lo, ymax=ci.95.hi),
								position=position_dodge(width=0.5),
								width=0.25, size=0.5) +
	geom_point(aes(shape=Policy), 
						 position=position_dodge(width=0.5),size=1.5) + 
	scale_y_continuous("% Support Policy",limits=c(0.2,0.83),
										 breaks=seq(0,1,0.2),
										 labels=100*seq(0,1,0.2)) + 
	scale_x_discrete("Subgroup",breaks=x_labels$Subgroup,labels=x_labels$Subgroup_n) + 
	facet_wrap("Policy",ncol = 1) + 
	# geom_text(data=desc_tab,aes(x = Subgroup,y=0.20,label=n)) +
	theme_bw() + 
	theme(legend.position = "none",
				strip.text=element_text(size=18),
				axis.text = element_text(size=12))

fname = "descriptive_points_panel_b.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_desc_points_b, width=11, height=5)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=11, height=5)
	print(plot_desc_points_b)
	dev.off()
}

fname = "descriptive_points_panel_b_withpers.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_desc_points_b, width=12.75, height=5)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=12.75, height=5)
	print(plot_desc_points_b)
	dev.off()
}

# panel c
x_labels <- desc_tab[which(desc_tab$variable=="dist_support_bi"),c("Subgroup","Subgroup_n")]
plot_desc_points_c <- ggplot(subset(desc_tab,variable=="dist_support_bi"),aes(x=Subgroup,y=mean)) + 
	geom_errorbar(aes(ymin=ci.95.lo, ymax=ci.95.hi),
								position=position_dodge(width=0.5),
								width=0.25, size=0.5) +
	geom_point(aes(shape=Policy), 
						 position=position_dodge(width=0.5),size=1.5) + 
	scale_y_continuous("% Support Policy",limits=c(0.2,0.83),
										 breaks=seq(0,1,0.2),
										 labels=100*seq(0,1,0.2)) + 
	scale_x_discrete("Subgroup",breaks=x_labels$Subgroup,labels=x_labels$Subgroup_n) + 
	facet_wrap("Policy",ncol = 1) + 
	# geom_text(data=desc_tab,aes(x = Subgroup,y=0.20,label=n)) +
	theme_bw() + 
	theme(legend.position = "none",
				strip.text=element_text(size=18),
				axis.text = element_text(size=12))

fname = "descriptive_points_panel_c.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_desc_points_c, width=11, height=5)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=11, height=5)
	print(plot_desc_points_c)
	dev.off()
}

fname = "descriptive_points_panel_c_withpers.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_desc_points_c, width=12.75, height=5)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=12.75, height=5)
	print(plot_desc_points_c)
	dev.off()
}

# format for table:
desc_tab_output <- bind_rows(desc_full,desc_pid,desc_ideo,desc_income,desc_death, desc_persexp)

desc_tab_output$Subgroup <- factor(desc_tab_output$Subgroup,
																	 levels = c("All\nRespondents",
																	 					 "Democrats","Independents","Republicans",
																	 					 # "Income <$50k","Income >$50k",
																	 					 "Liberals","Moderates","Conservatives",
																	 					 
																	 					 "Below Median\nIncome",
																	 					 "Above Median\nIncome",
																	 					 "Below Median\nOverdose Rate",
																	 					 "Above Median\nOverdose Rate",
																	 					 "Know\nSomeone\nw/ Addiction","Don't Know\nSomeone\nw/ Addiction"#,
																	 ),
																	 ordered = T)
desc_tab_output$Policy <- factor(desc_tab_output$variable,levels = c("ability_support_bi","needs_support_bi","dist_support_bi"),
																 labels=c("Treatment Funding:\nRedistributive","Treatment Funding:\nNeeds-Based","Clinic Construction"),
																 ordered=T)
desc_tab_output$sd_text <- paste0("(",round(desc_tab_output$sd,3),")")

desc_tab_output <- desc_tab_output %>%
	select(Subgroup,Policy,mean,sd_text,n) %>%
	gather(value="value",key="stat",-Policy,-Subgroup)

desc_tab_output1 <- desc_tab_output %>%
	filter(Policy=="Treatment Funding:\nRedistributive") %>%
	spread(key="stat",value="value",convert = T) %>%
	select(Subgroup,mean,sd_text,n)
desc_tab_output2 <- desc_tab_output %>%
	filter(Policy=="Treatment Funding:\nNeeds-Based") %>%
	spread(key="stat",value="value",convert = T) %>%
	select(mean,sd_text,n)
desc_tab_output3 <- desc_tab_output %>%
	filter(Policy=="Clinic Construction") %>%
	spread(key="stat",value="value",convert = T) %>%
	select(mean,sd_text,n)

desc_tab_output <- bind_cols(desc_tab_output1,desc_tab_output2,desc_tab_output3)
colnames(desc_tab_output) <- c("Subgroup","Mean","SD","n","Mean","SD","n","Mean","SD","n")

## Table A4:
print(xtable(desc_tab_output,digits = 3),
			include.rownames=F,
			add.to.row = list(pos=list(-1),command=" & & Redistributive & & & Overdose rate-based & & & Clinic Construction & \\\\ "),
			booktabs=T,file="means_tab.tex",floating = F
)

nimby_pid <- lm(dist_support_bi ~ factor(pid) ,data=a);summary(nimby_pid)
nimby_pid_controls <- lm(dist_support_bi ~ factor(pid) + AGE + INCOME + female + dem + rep + liberal + conservative + college + married + white + black + hispanic + own,data=a);summary(nimby_pid_controls)
writeLines(as.character(100*abs(round(nimby_pid$coefficients[3],2))), "Tables/nimby_pid_perc.tex",sep = "")


# --------------------------------- #
#### Analyses: Treatment Funding ####
# --------------------------------- #

fit_redist <- lm(ability_support_bi ~ above_income + above_death_rate + AGE + female + black + hispanic + own  + factor(pid) + factor(ideo) + college + married , data=a)
fit_needs <- lm(needs_support_bi ~ above_income + above_death_rate + AGE + female + black + hispanic + own  + factor(pid) + factor(ideo) + college + married, data=a)
fit_redist_nocov <- lm(ability_support_bi ~ above_income, data=a)
fit_needs_nocov <- lm(needs_support_bi ~ above_death_rate, data=a)

fit_redist_pid_nocov <- lm(ability_support_bi ~ factor(pid), data=a)
fit_needs_pid_nocov <- lm(needs_support_bi ~ factor(pid), data=a)

fit_redist_intx_nocov <- lm(ability_support_bi ~ above_income*factor(pid), data=a)
fit_needs_intx_nocov <- lm(needs_support_bi ~ above_death_rate*factor(pid), data=a)

fit_redist_rep_nocov <- lm(ability_support_bi ~ above_income, data=subset(a,pid==1))
fit_needs_rep_nocov <- lm(needs_support_bi ~ above_death_rate, data=subset(a,pid==1))
fit_redist_dem_nocov <- lm(ability_support_bi ~ above_income, data=subset(a,pid==-1))
fit_needs_dem_nocov <- lm(needs_support_bi ~ above_death_rate, data=subset(a,pid==-1))

fit_redist_rep <- lm(ability_support_bi ~ AGE + female + black + hispanic + own  + factor(ideo) + college + married + above_income + above_death_rate, data=subset(a,pid==1))
fit_needs_rep <- lm(needs_support_bi ~ AGE + female + black + hispanic + own + factor(ideo) + college + married + above_income + above_death_rate, data=subset(a,pid==1))
fit_redist_dem <- lm(ability_support_bi ~ AGE + female + black + hispanic + own + factor(ideo) + college + married + above_income + above_death_rate, data=subset(a,pid==-1))
fit_needs_dem <- lm(needs_support_bi ~ AGE + female + black + hispanic + own + factor(ideo) + college + married + above_income + above_death_rate, data=subset(a,pid==-1))



#### Output effects in percentage point terms: 
writeLines(as.character(abs(100*round(fit_redist_nocov$coefficients["above_income"],2))), "selfinterest_effect_redist_perc.tex",sep = "")
writeLines(as.character(abs(100*round(fit_needs_nocov$coefficients["above_death_rate"],2))), "selfinterest_effect_needs_perc.tex",sep = "")
writeLines(as.character(abs(100*round(fit_redist_pid_nocov$coefficients["factor(pid)1"],2))), "rep_effect_redist_perc.tex",sep = "")
writeLines(as.character(abs(100*round(fit_needs_pid_nocov$coefficients["factor(pid)1"],2))), "rep_effect_needs_perc.tex",sep = "")

## differential effects:
writeLines(as.character(abs(100*round(fit_redist_intx_nocov$coefficients["above_income"] + fit_redist_intx_nocov$coefficients["above_income:factor(pid)1"],2))), "selfinterest_effect_rep_redist_perc.tex",sep = "")

## interactive effects:
writeLines(as.character(abs(100*round(fit_redist_intx_nocov$coefficients["above_income:factor(pid)1"],2))), "selfinterest_intx_redist_perc.tex",sep = "")


## Table A5:
stargazer(fit_redist_nocov,fit_redist,fit_redist_pid_nocov,fit_needs_nocov,fit_needs,fit_needs_pid_nocov,
					dep.var.caption = "Support for:",
					dep.var.labels = c("Redistributive Policy","Needs-based Policy"),
					covariate.labels = c("Above State Median Income",
															 "Above Median Overdose Rate",
															 "Age",
															 "Female",
															 "Black",
															 "Hispanic",
															 "Homeowner",
															 "Independent",
															 "Republican",
															 "Moderate",
															 "Conservative",
															 "College Degree",
															 "Married"),
					notes = c("Omitted category for partisanship is 'Democrat'","and for ideology is 'Liberal'"),
					notes.align = "l",
					float = F,df=F,omit.stat = c("rsq","ser"),out = "reg_tab_burden.tex"
)

## interactive effects:

redist_intx_rep <- lm(ability_support_bi ~ above_income*factor(pid),data=a);summary(redist_intx_rep)
redist_intx_rep_controls <- lm(ability_support_bi ~ above_income*factor(pid) + AGE + INCOME + female + liberal + conservative + college + married + white + black + hispanic + own + above_death_rate, data=a); summary(redist_intx_rep_controls)

redist_intx_death <- lm(ability_support_bi ~ above_income*above_death_rate,data=a);summary(redist_intx_death)
redist_intx_death_controls <- lm(ability_support_bi ~ above_income*above_death_rate + AGE + INCOME + female + dem + rep + liberal + conservative + college + married + white + black + hispanic + own, data=a); summary(redist_intx_death_controls)

needs_intx_rep <- lm(needs_support_bi ~ above_death_rate*factor(pid),data=a);summary(needs_intx_rep)
needs_intx_rep_controls <- lm(needs_support_bi ~ above_death_rate*factor(pid) + AGE + INCOME + female + liberal + conservative + college + married + white + black + hispanic + own + above_income, data=a); summary(needs_intx_rep_controls)

needs_intx_inc <- lm(needs_support_bi ~ above_income*above_death_rate,data=a);summary(needs_intx_inc)
needs_intx_inc_controls <- lm(needs_support_bi ~ above_income*above_death_rate + AGE + INCOME + female + dem + rep + liberal + conservative + college + married + white + black + hispanic + own, data=a); summary(needs_intx_inc_controls)


## Table A6: heterogeneous treatment effects
stargazer(
	# redist_intx_death,redist_intx_death_controls,
	redist_intx_rep,redist_intx_rep_controls,
	# needs_intx_inc,needs_intx_inc_controls,
	needs_intx_rep,needs_intx_rep_controls,
	omit = c("AGE", "INCOME", "female", "dem", "rep", "liberal", "conservative", "college", "married", "white", "black", "hispanic", "own"),
	covariate.labels = c(
		"Above Median Income", 
		"Party = Independent",
		"Party = Republican",
		"Above Median Overdose Rate $\\times$ Party = Independent",
		"Above Median Overdose Rate $\\times$ Party = Republican",
		"Above Median Overdose Rate",
		"Above Median Income $\\times$ Party = Independent",
		"Above Median Income $\\times$ Party = Republican"),
	dep.var.labels = c("Redistributive Policy","Needs-based Policy"),dep.var.caption = "Support for:",
	omit.labels = rep("Demographic controls",13), # have to go in and manually delete extra rows added to table here
	omit.yes.no = c("$\\checkmark$",""),
	# add.lines = list(
	# 	c("Demographic controls", "", "$\\checkmark$", "", "$\\checkmark$ ")
	# ),
	notes = "Omitted category for partisanship is 'Democrat'",
	float = F,df=F,omit.stat = c("rsq","ser"),out = "treat_burden_tab_controls_v2.tex"
)


# generate master plot of party ID effects for both treatment funding proposals (for Figure A2)
treatsumm_pid <- subset(a, !is.na(support_bi) & !is.na(pid) & pid!=0) %>%
	group_by(pid,condition,above_death_rate, above_income) %>%
	summarize(groupmean = mean(support_bi),
						groupsd = sd(support_bi),
						n = n()) %>%
	mutate(ci.95.hi = groupmean + qnorm(0.975)*groupsd/sqrt(n),
				 ci.95.lo = groupmean + qnorm(0.025)*groupsd/sqrt(n)) 

treatsumm_pid$type<-ifelse(treatsumm_pid$above_death_rate==1 & treatsumm_pid$above_income==1, "High Income,\nHigh Overdose Rate", 
													 ifelse(treatsumm_pid$above_death_rate==0 & treatsumm_pid$above_income==1, "High Income,\nLow Overdose Rate",
													 			 ifelse(treatsumm_pid$above_death_rate==1 & treatsumm_pid$above_income==0, "Low Income,\nHigh Overdose Rate",
													 			 			 "Low Income,\nLow Overdose Rate")))

treatsumm_pid_redist <- subset(a, !is.na(support_bi) & !is.na(pid) & pid!=0) %>%
	group_by(pid,condition, above_income) %>%
	summarize(groupmean = mean(support_bi),
						groupsd = sd(support_bi),
						n = n()) %>%
	mutate(ci.95.hi = groupmean + qnorm(0.975)*groupsd/sqrt(n),
				 ci.95.lo = groupmean + qnorm(0.025)*groupsd/sqrt(n)) 

treatsumm_pid_needs <- subset(a, !is.na(support_bi) & !is.na(pid) & pid!=0) %>%
	group_by(pid,condition, above_death_rate) %>%
	summarize(groupmean = mean(support_bi),
						groupsd = sd(support_bi),
						n = n()) %>%
	mutate(ci.95.hi = groupmean + qnorm(0.975)*groupsd/sqrt(n),
				 ci.95.lo = groupmean + qnorm(0.025)*groupsd/sqrt(n)) 

treatsumm_pid_both <- bind_rows(treatsumm_pid_redist,treatsumm_pid_needs)

treatsumm_pid_both$type <- ifelse(!is.na(treatsumm_pid_both$above_death_rate),ifelse(treatsumm_pid_both$above_death_rate==1, "High Overdose Rate", "Low Overdose Rate"),
																	ifelse(treatsumm_pid_both$above_income==0, "Low Income","High Income"))

# treatsumm_pid_both

## Figure A2:
plot_pid_points <- ggplot(subset(treatsumm_pid_both,(condition=="ability" & !is.na(above_income)) | (condition=="need" & !is.na(above_death_rate))),aes(group=factor(pid))) + 
	geom_errorbar(aes(ymin=ci.95.lo, ymax=ci.95.hi,x=type,color=factor(pid)),
								position=position_dodge(width=0.5),
								width=.1, size=0.5) +
	geom_point(aes(x=type,y=groupmean,color=factor(pid),shape=factor(pid)),size=2,
						 position=position_dodge(width=0.5)) + 
	facet_wrap("condition",scales = "free_x",
						 labeller = labeller(condition = c(`ability` = "Income-Based", `need` = "Overdose Rate-Based"))
	) + 
	scale_y_continuous("% Support Policy",limits=c(0.15,0.85),
										 breaks=seq(0,1,0.2),
										 labels=100*seq(0,1,0.2)) + 
	scale_color_manual(breaks=c(-1,1),labels=c("Democrats","Republicans"),values=c("blue","red")) + 
	scale_shape_manual(breaks=c(-1,1),labels=c("Democrats","Republicans"),values=c(16,17)) + 
	# scale_shape_manual(breaks=c("ability", "need"),labels=c("Ability","Need"),values=c(15,17)) + 
	scale_x_discrete("Subgroup") +
	theme_bw() + 
	theme(legend.position = c(0.9,0.82)) +
	theme(legend.title=element_blank())

## Output:
fname = "burden_subgroup_points.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_pid_points, width=7, height=7)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=8, height=5)
	print(plot_pid_points)
	dev.off()
}

### Effect of financial self-interest:
treatment_tab <- data.frame(
	rbind(
		c(summary(fit_redist_nocov)$coefficients["above_income",c("Estimate","Std. Error")],n=nrow(model.frame(fit_redist_nocov))),
		c(summary(fit_redist)$coefficients["above_income",c("Estimate","Std. Error")],n=nrow(model.frame(fit_redist))),
		c(summary(fit_needs_nocov)$coefficients["above_death_rate",c("Estimate","Std. Error")],n=nrow(model.frame(fit_needs_nocov))),
		c(summary(fit_needs)$coefficients["above_death_rate",c("Estimate","Std. Error")],n=nrow(model.frame(fit_needs)))
	)
)
treatment_tab <- mutate(treatment_tab,
												ci.lo.95 = Estimate + qnorm(0.025)*Std..Error,
												ci.hi.95 = Estimate + qnorm(0.975)*Std..Error,
												ci.lo.90 = Estimate + qnorm(0.05)*Std..Error,
												ci.hi.90 = Estimate + qnorm(0.95)*Std..Error
)

treatment_tab$cov <- factor(c(0,1,0,1))
treatment_tab$frame <- c("Redistributive","Redistributive","Overdose\nrate-based","Overdose\nrate-based")
treatment_tab$frame <- factor(treatment_tab$frame,levels = c("Redistributive","Overdose\nrate-based"),ordered = T)

treatment_tab

## Figure 2:
plot_treat_selfinterest <- ggplot(treatment_tab, aes(x=frame,group=cov)) +
	geom_errorbar(aes(ymin=ci.lo.95, ymax=ci.hi.95),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=0.35) +
	geom_errorbar(aes(ymin=ci.lo.90, ymax=ci.hi.90),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=1.5) +
	geom_point(aes(y=Estimate,shape=cov),position=position_dodge(width=0.2),
						 data=treatment_tab,size=4) + 
	geom_hline(aes(yintercept=0),lty=2) + 
	# xlab(NULL) +
	ylab("Effect of financial self-interest\non policy support") +
	theme_bw() + 
	# scale_y_continuous(limits=c(-0.175, 0.005)) + 
	scale_x_discrete("Funding model") + 
	scale_shape_discrete("Covariate-\nadjusted",labels=c("No","Yes")) + 
	# coord_flip() + 
	theme(text=element_text(colour=black,
													# family="CM Roman",
													size=15),
				legend.position = "right")

## Output:
fname = "selfinterest_effects_burden.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_treat_selfinterest, width=6, height=4)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=6, height=4)
	print(plot_treat_selfinterest)
	dev.off()
}

## Effects of PID:
treatment_tab <- data.frame(
	rbind(
		c(summary(fit_redist_pid_nocov)$coefficients["factor(pid)1",c("Estimate","Std. Error")],n=nrow(model.frame(fit_redist_pid_nocov))),
		c(summary(fit_redist)$coefficients["factor(pid)1",c("Estimate","Std. Error")],n=nrow(model.frame(fit_redist))),
		c(summary(fit_needs_pid_nocov)$coefficients["factor(pid)1",c("Estimate","Std. Error")],n=nrow(model.frame(fit_needs_pid_nocov))),
		c(summary(fit_needs)$coefficients["factor(pid)1",c("Estimate","Std. Error")],n=nrow(model.frame(fit_needs)))
	)
)
treatment_tab <- mutate(treatment_tab,
												ci.lo.95 = Estimate + qnorm(0.025)*Std..Error,
												ci.hi.95 = Estimate + qnorm(0.975)*Std..Error,
												ci.lo.90 = Estimate + qnorm(0.05)*Std..Error,
												ci.hi.90 = Estimate + qnorm(0.95)*Std..Error
)

treatment_tab$cov <- factor(c(0,1,0,1))
treatment_tab$frame <- c("Redistributive","Redistributive","Overdose\nrate-based","Overdose\nrate-based")
treatment_tab$frame <- factor(treatment_tab$frame,levels = c("Redistributive","Overdose\nrate-based"),ordered = T)

treatment_tab

## Figure 3:
plot_treat_rep <- ggplot(treatment_tab, aes(x=frame,group=cov)) +
	geom_errorbar(aes(ymin=ci.lo.95, ymax=ci.hi.95),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=0.35) +
	geom_errorbar(aes(ymin=ci.lo.90, ymax=ci.hi.90),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=1.5) +
	geom_point(aes(y=Estimate,shape=cov),position=position_dodge(width=0.2),
						 data=treatment_tab,size=4) + 
	geom_hline(aes(yintercept=0),lty=2) + 
	# xlab(NULL) +
	ylab("Effect of partisanship\non policy support") +
	theme_bw() + 
	# scale_y_continuous(limits=c(-0.175, 0.005)) + 
	scale_x_discrete("Funding model") + 
	scale_shape_discrete("Covariate-\nadjusted",labels=c("No","Yes")) + 
	# coord_flip() + 
	theme(text=element_text(colour=black,
													# family="CM Roman",
													size=15),
				legend.position = "right")

# Output:
fname = "rep_effects_burden.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_treat_rep, width=6, height=4)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=6, height=4)
	print(plot_treat_rep)
	dev.off()
}

## Effect of financial self-interest interacted with PID:
treatment_tab <- data.frame(
	rbind(
		c(summary(fit_redist_rep_nocov)$coefficients["above_income",c("Estimate","Std. Error")],n=nrow(model.frame(fit_redist_rep_nocov))),
		c(summary(fit_redist_rep)$coefficients["above_income",c("Estimate","Std. Error")],n=nrow(model.frame(fit_redist_rep))),
		c(summary(fit_needs_rep_nocov)$coefficients["above_death_rate",c("Estimate","Std. Error")],n=nrow(model.frame(fit_needs_rep_nocov))),
		c(summary(fit_needs_rep)$coefficients["above_death_rate",c("Estimate","Std. Error")],n=nrow(model.frame(fit_needs_rep))),
		c(summary(fit_redist_dem_nocov)$coefficients["above_income",c("Estimate","Std. Error")],n=nrow(model.frame(fit_redist_dem_nocov))),
		c(summary(fit_redist_dem)$coefficients["above_income",c("Estimate","Std. Error")],n=nrow(model.frame(fit_redist_dem))),
		c(summary(fit_needs_dem_nocov)$coefficients["above_death_rate",c("Estimate","Std. Error")],n=nrow(model.frame(fit_needs_dem_nocov))),
		c(summary(fit_needs_dem)$coefficients["above_death_rate",c("Estimate","Std. Error")],n=nrow(model.frame(fit_needs_dem)))
	)
)
treatment_tab <- mutate(treatment_tab,
												ci.lo.95 = Estimate + qnorm(0.025)*Std..Error,
												ci.hi.95 = Estimate + qnorm(0.975)*Std..Error,
												ci.lo.90 = Estimate + qnorm(0.05)*Std..Error,
												ci.hi.90 = Estimate + qnorm(0.95)*Std..Error
)

treatment_tab$cov <- factor(c(0,1,0,1,0,1,0,1))
treatment_tab$Partisanship <- c("Republicans","Republicans","Republicans","Republicans","Democrats","Democrats","Democrats","Democrats")
treatment_tab$frame <- c("Redistributive","Redistributive","Overdose\nrate-based","Overdose\nrate-based","Redistributive","Redistributive","Overdose\nrate-based","Overdose\nrate-based")
treatment_tab$frame <- factor(treatment_tab$frame,levels = c("Redistributive","Overdose\nrate-based"),ordered = T)

treatment_tab

## Figure 4:
plot_treat_selfinterest_intx <- ggplot(treatment_tab, aes(x=Partisanship,group=cov)) +
	geom_errorbar(aes(ymin=ci.lo.95, ymax=ci.hi.95),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=0.35) +
	geom_errorbar(aes(ymin=ci.lo.90, ymax=ci.hi.90),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=1.5) +
	geom_point(aes(y=Estimate,shape=cov),position=position_dodge(width=0.2),
						 data=treatment_tab,size=4) + 
	geom_hline(aes(yintercept=0),lty=2) + 
	facet_wrap("frame") + 
	# xlab(NULL) +
	ylab("Effect of financial self-interest\non policy support") +
	theme_bw() + 
	# scale_y_continuous(limits=c(-0.175, 0.005)) + 
	scale_x_discrete("Partisanship",limits=c("Republicans","Democrats")) + 
	scale_shape_discrete("Covariate-\nadjusted",labels=c("No","Yes")) + 
	# coord_flip() + 
	theme(text=element_text(colour=black,
													# family="CM Roman",
													size=15),
				legend.position = "right")

fname = "selfinterest_intx_burden.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_treat_selfinterest_intx, width=9, height=6)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=9, height=6)
	print(plot_treat_selfinterest_intx)
	dev.off()
}


#-----------------------#
#### Analyses: NIMBY ####
#-----------------------#
# Experimental treatment effect:
fit_dist_nocov <- lm(dist_support_bi ~ distance, data=a)
fit_dist <- lm(dist_support_bi ~ distance + above_income + above_death_rate + AGE + female + black + hispanic + own  + factor(pid) + factor(ideo) + college + married, data=a)

## Table A7:
stargazer(fit_dist_nocov,fit_dist,
					dep.var.caption = "Support for:",
					dep.var.labels = c("Clinic Construction"),
					covariate.labels = c("Distance Condition",
															 "Above State Median Income",
															 "Above Median Overdose Rate",
															 "Age",
															 "Female",
															 "Black",
															 "Hispanic",
															 "Homeowner",
															 "Independent",
															 "Republican",
															 "Moderate",
															 "Conservative",
															 "College Degree",
															 "Married"),
					notes = c("Omitted category for partisanship is 'Democrat'","and for ideology is 'Liberal'"),
					notes.align = "l",
					float = F,df=F,omit.stat = c("rsq","ser"),out = "reg_tab_nimby.tex"
)



# create factor labels for plotting:
a$distance <- factor(a$distance,
										 levels=c("far","near"),
										 labels=c("Distance:\nFar","Distance:\nNear"))


distance.90 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear"],a$dist_support_bi[a$distance=="Distance:\nFar"],conf.level = 0.9)
distance.95 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear"],a$dist_support_bi[a$distance=="Distance:\nFar"],conf.level = 0.95)


# output treatment effect:
writeLines(as.character(round(distance.95$estimate[2]-distance.95$estimate[1],2)*100),"nimby_treatment_effect_perc.tex",sep = "")
writeLines(as.character(round(distance.95$estimate[2],2)*100),"nimby_support_far_perc.tex",sep = "")
writeLines(as.character(round(distance.95$estimate[1],2)*100),"nimby_support_near_perc.tex",sep = "")


# treatment effects in subgroups with and without controls:

nimby_full_nocov <- lm(dist_support_bi ~ distance ,data=subset(a))
nimby_rep_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,pid==1))
nimby_dem_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,pid==-1))
nimby_cons_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,ideo==1))
nimby_lib_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,ideo==-1))
nimby_white_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,white==1))
nimby_black_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,black==1))
nimby_hiinc_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,above_income==1))
nimby_loinc_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,above_income==0))
nimby_hineed_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,above_death_rate==1))
nimby_loneed_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,above_death_rate==0))
nimby_pers1_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,personal_bi==1))
nimby_pers0_nocov <- lm(dist_support_bi ~ distance ,data=subset(a,personal_bi==0))

nimby_full <- lm(dist_support_bi ~ distance + AGE + female + black + hispanic + own + factor(pid) + factor(ideo) + college + married + above_income + above_death_rate, data=subset(a))
nimby_rep <- lm(dist_support_bi ~ distance  + AGE + female + black + hispanic + own + factor(ideo) + college + married + above_income + above_death_rate,data=subset(a,pid==1))
nimby_dem <- lm(dist_support_bi ~ distance  + AGE + female + black + hispanic + own + factor(ideo) + college + married + above_income + above_death_rate,data=subset(a,pid==-1))
nimby_cons <- lm(dist_support_bi ~ distance  + AGE + female + black + hispanic + own + factor(pid) + college + married + above_income + above_death_rate,data=subset(a,ideo==1))
nimby_lib <- lm(dist_support_bi ~ distance  + AGE + female + black + hispanic + own + factor(pid) + college + married + above_income + above_death_rate,data=subset(a,ideo==-1))
nimby_white <- lm(dist_support_bi ~ distance  + AGE + female + own + factor(pid) + factor(ideo) + college + married + above_income + above_death_rate,data=subset(a,white==1))
nimby_black <- lm(dist_support_bi ~ distance  + AGE + female + own + factor(pid) + factor(ideo) + college + married + above_income + above_death_rate,data=subset(a,black==1))
nimby_hiinc <- lm(dist_support_bi ~ distance + AGE + female + black + hispanic + own + factor(pid) + factor(ideo) + college + married +  above_death_rate, data=subset(a,above_income==1))
nimby_loinc <- lm(dist_support_bi ~ distance + AGE + female + black + hispanic + own + factor(pid) + factor(ideo) + college + married + above_death_rate, data=subset(a,above_income==0))
nimby_hineed <- lm(dist_support_bi ~ distance + AGE + female + black + hispanic + own + factor(pid) + factor(ideo) + college + married + above_income, data=subset(a,above_death_rate==1))
nimby_loneed <- lm(dist_support_bi ~ distance + AGE + female + black + hispanic + own + factor(pid) + factor(ideo) + college + married + above_income, data=subset(a,above_death_rate==0))
nimby_pers1 <- lm(dist_support_bi ~ distance + AGE + female + black + hispanic + own + factor(pid) + factor(ideo) + college + married + above_income + above_death_rate, data=subset(a,personal_bi==1))
nimby_pers0 <- lm(dist_support_bi ~ distance + AGE + female + black + hispanic + own + factor(pid) + factor(ideo) + college + married + above_income + above_death_rate, data=subset(a,personal_bi==0))


nimby_intx_pid <- lm(dist_support_bi ~ distance*factor(pid),data=a)
nimby_intx_pid_controls <- lm(dist_support_bi ~ distance*factor(pid) + AGE + INCOME + female + liberal + conservative + college + married + white + black + hispanic + own,data=a)


nimby_intx_inc <- lm(dist_support_bi ~ distance*above_income,data=a)
nimby_intx_inc_controls <- lm(dist_support_bi ~ distance*above_income + AGE + INCOME + female + dem + rep + liberal + conservative + college + married + white + black + hispanic + own,data=a)

nimby_intx_need <- lm(dist_support_bi ~ distance*above_death_rate,data=a)
nimby_intx_need_controls <- lm(dist_support_bi ~ distance*above_death_rate + AGE + INCOME + female + dem + rep + liberal + conservative + college + married + white + black + hispanic + own,data=a)



# based on personal exposure
distance_pers1.90 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear" & a$personal_bi==1],a$dist_support_bi[a$distance=="Distance:\nFar" & a$personal_bi==1],conf.level = 0.9)
distance_pers1.95 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear" & a$personal_bi==1],a$dist_support_bi[a$distance=="Distance:\nFar" & a$personal_bi==1],conf.level = 0.95)
distance_pers0.90 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear" & a$personal_bi==0],a$dist_support_bi[a$distance=="Distance:\nFar" & a$personal_bi==0],conf.level = 0.9)
distance_pers0.95 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear" & a$personal_bi==0],a$dist_support_bi[a$distance=="Distance:\nFar" & a$personal_bi==0],conf.level = 0.95)

nimby_intx_pers <- lm(dist_support_bi ~ distance*personal_bi,data=a);summary(nimby_intx_pers) # neg insig intx
nimby_intx_pers_controls <- lm(dist_support_bi ~ distance*personal_bi + AGE + INCOME + female + dem + rep + liberal + conservative + college + married + white + black + hispanic + own,data=a);summary(nimby_intx_pers_controls) # neg insig intx



## based on population density:
distance_dens1.90 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear" & a$pop_density_bin_num>2],a$dist_support_bi[a$distance=="Distance:\nFar" & a$pop_density_bin_num>2],conf.level = 0.9)
distance_dens1.95 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear" & a$pop_density_bin_num>2],a$dist_support_bi[a$distance=="Distance:\nFar" & a$pop_density_bin_num>2],conf.level = 0.95)
distance_dens0.90 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear" & a$pop_density_bin_num<3],a$dist_support_bi[a$distance=="Distance:\nFar" & a$pop_density_bin_num<3],conf.level = 0.9)
distance_dens0.95 <- t.test(a$dist_support_bi[a$distance=="Distance:\nNear" & a$pop_density_bin_num<3],a$dist_support_bi[a$distance=="Distance:\nFar" & a$pop_density_bin_num<3],conf.level = 0.95)

nimby_intx_density <- lm(dist_support_bi ~ distance*pop_density,data=a) # small sig pos intx
nimby_intx_density_bin <- lm(dist_support_bi ~ distance*pop_density_bin_num,data=a) # insig pos intx
nimby_intx_density_01 <- lm(dist_support_bi ~ distance*(pop_density_bin_num>2),data=a) # insig neg intx



## Table A9:
stargazer(nimby_intx_inc,nimby_intx_inc_controls,
					nimby_intx_need,nimby_intx_need_controls,
					nimby_intx_pers,nimby_intx_pers_controls,
					nimby_intx_pid,nimby_intx_pid_controls,
					omit = c("AGE", "INCOME", "female", "dem", "rep", "liberal", "conservative", "college", "married", "white", "black", "hispanic", "own"),
					covariate.labels = c("Treatment Condition = Near", 
															 "Above Median Income", 
															 "Treatment Condition = Near $\\times$ Above Median Income",
															 "Above Median Overdose Rate",
															 "Treatment Condition = Near $\\times$ Above Median Overdose Rate",
															 "Know Someone w/ Addiction",
															 "Treatment Condition = Near $\\times$ Know Someone w/ Addiction",
															 "Party = Independent", "Party = Republican",
															 "Treatment Condition = Near $\\times$ Party = Independent",
															 "Treatment Condition = Near $\\times$ Party = Republican"),
					dep.var.labels = "Clinic Construction Support",
					omit.labels = rep("Demographic controls",13), # have to manually delete extra lines added here
					omit.yes.no = c("$\\checkmark$",""),
					notes = c("Omitted category for partisanship is 'Democrat'. Controls include age, income, gender,","partisanship, ideology, college degree, marital status, race, and homeownership."),notes.align = "l",
					float = F,df=F,omit.stat = c("rsq","ser"),out = "treat_nimby_tab_controls.tex"
)



# summary of mean responses by respondent characteristics:
treatsumm_inc <- subset(a, !is.na(dist_support_bi) & !is.na(distance)) %>%
	group_by(distance,above_income) %>%
	summarize(groupmean = mean(dist_support_bi),
						groupsd = sd(dist_support_bi),
						n = n()) %>%
	mutate(ci.95.hi = groupmean + qnorm(0.975)*groupsd/sqrt(n),
				 ci.95.lo = groupmean + qnorm(0.025)*groupsd/sqrt(n)) 

## Figure A4:
plot_inc_points <- ggplot(treatsumm_inc,aes(group=distance)) + 
	geom_errorbar(aes(ymin=ci.95.lo, ymax=ci.95.hi,x=above_income),
								position=position_dodge(width=0.5),
								width=.1, size=0.5) +
	geom_point(data=treatsumm_inc,aes(x=above_income,y=groupmean,shape=distance),size=1.5,
						 position=position_dodge(width=0.5)) + 
	scale_y_continuous("% Support Policy",limits=c(0.2,0.8),
										 breaks=seq(0,1,0.2),
										 labels=100*seq(0,1,0.2)) + 
	scale_x_continuous("Income",breaks=c(0,1),labels=c("Below average","Above average")) + 
	scale_shape_discrete("Distance Treatment",breaks=c("Distance:\nFar","Distance:\nNear"),labels=c("Far","Near")) +
	theme_bw() + 
	theme(legend.position = c(0.7,0.85))

## Output:
fname = "nimby_inc_effect_points.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_inc_points, width=5, height=4)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=5, height=4)
	print(plot_inc_points)
	dev.off()
}


# state death rate 
treatsumm_need <- subset(a, !is.na(dist_support_bi) & !is.na(distance)) %>%
	group_by(distance,above_death_rate) %>%
	summarize(groupmean = mean(dist_support_bi),
						groupsd = sd(dist_support_bi),
						n = n()) %>%
	mutate(ci.95.hi = groupmean + qnorm(0.975)*groupsd/sqrt(n),
				 ci.95.lo = groupmean + qnorm(0.025)*groupsd/sqrt(n)) 

## Figure A5:
plot_need_points <- ggplot(treatsumm_need,aes(group=distance)) + 
	geom_errorbar(aes(ymin=ci.95.lo, ymax=ci.95.hi,x=above_death_rate),
								position=position_dodge(width=0.5),
								width=.1, size=0.5) +
	geom_point(data=treatsumm_need,aes(x=above_death_rate,y=groupmean,shape=distance),size=1.5,
						 position=position_dodge(width=0.5)) + 
	scale_y_continuous("% Support Policy",limits=c(0.2,0.8),
										 breaks=seq(0,1,0.2),
										 labels=100*seq(0,1,0.2)) + 
	scale_x_continuous("Overdose Rate",breaks=c(0,1),labels=c("Below median","Above median")) + 
	scale_shape_discrete("Distance Treatment",breaks=c("Distance:\nFar","Distance:\nNear"),labels=c("Far","Near")) +
	theme_bw() + 
	theme(legend.position = c(0.7,0.85))

## Output:
fname = "nimby_need_effect_points.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_need_points, width=7, height=7)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=5, height=4)
	print(plot_need_points)
	dev.off()
}

treatsumm_pid <- subset(a, !is.na(dist_support_bi) & !is.na(distance) & !is.na(pid)) %>%
	group_by(distance,pid) %>%
	summarize(groupmean = mean(dist_support_bi),
						groupsd = sd(dist_support_bi),
						n = n()) %>%
	mutate(ci.95.hi = groupmean + qnorm(0.975)*groupsd/sqrt(n),
				 ci.95.lo = groupmean + qnorm(0.025)*groupsd/sqrt(n)) 

## Figure A6:
plot_pid_points <- ggplot(treatsumm_pid) + 
	geom_errorbar(aes(ymin=ci.95.lo, ymax=ci.95.hi,x=distance,group=factor(pid),color=factor(pid)),
								position=position_dodge(width=0.5),
								width=.1, size=0.5) +
	geom_point(data=treatsumm_pid,aes(x=distance,y=groupmean,shape=factor(pid),color=factor(pid)),size=1.5,
						 position=position_dodge(width=0.5)) + 
	scale_y_continuous("% Support Policy",limits=c(0.15,0.85),
										 breaks=seq(0,1,0.2),
										 labels=100*seq(0,1,0.2)) + 
	scale_color_manual("Party ID",breaks=c(-1,0,1),labels=c("Democrats","Independents","Republicans"),values=c("blue","black","red")) + 
	scale_shape_manual("Party ID",breaks=c(-1,0,1),labels=c("Democrats","Independents","Republicans"),values=c(15,16,17)) + 
	scale_x_discrete("Distance",breaks=c("Distance:\nFar","Distance:\nNear"),labels=c("Far","Near")) +
	theme_bw() + 
	theme(legend.position = c(0.8,0.82))

## Output:
fname = "nimby_pid_effect_points.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_pid_points, width=7, height=7)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=5, height=4)
	print(plot_pid_points)
	dev.off()
}

# by personal exposure
summ_pers <- subset(a, !is.na(dist_support_bi) & !is.na(distance)) %>%
	group_by(distance,personal_bi) %>%
	summarize(groupmean = mean(dist_support_bi),
						groupsd = sd(dist_support_bi),
						n = n()) %>%
	mutate(ci.95.hi = groupmean + qnorm(0.975)*groupsd/sqrt(n),
				 ci.95.lo = groupmean + qnorm(0.025)*groupsd/sqrt(n))  %>%
	na.omit()

# summ_pers

## Figure A7:
plot_nimby_points_pers <- ggplot(summ_pers,aes(group=distance)) + 
	geom_errorbar(aes(ymin=ci.95.lo, ymax=ci.95.hi,x=personal_bi),
								position=position_dodge(width=0.5),
								width=.1, size=0.5) +
	geom_point(data=summ_pers,aes(x=personal_bi,y=groupmean,shape=distance),size=1.5,
						 position=position_dodge(width=0.5)) + 
	scale_y_continuous("% Support Policy",limits=c(0.2,0.8),
										 breaks=seq(0,1,0.2),
										 labels=100*seq(0,1,0.2)) + 
	scale_x_discrete("Know Someone w/ Addiction",limits=c(0,1),labels=c("No","Yes")) + 
	scale_shape_discrete("Distance Treatment",breaks=c("Distance:\nFar","Distance:\nNear"),labels=c("Far","Near")) +
	theme_bw() + 
	theme(legend.position = c(0.7,0.85))

## Output:
fname = "nimby_points_personal.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_nimby_points_pers, width=7, height=7)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=5, height=4)
	print(plot_nimby_points_pers)
	dev.off()
}
summary(lm(dist_support_bi ~ distance*personal_bi, a)) # insig intx

## Treatment effects for full sample and partisan subgroups:
treatment_tab <- data.frame(
	rbind(
		c(summary(nimby_full_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_full_nocov))),
		c(summary(nimby_full)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_full))),
		c(summary(nimby_rep_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_rep_nocov))),
		c(summary(nimby_rep)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_rep))),
		c(summary(nimby_dem_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_dem_nocov))),
		c(summary(nimby_dem)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_dem)))
	)
)
treatment_tab <- mutate(treatment_tab,
												ci.lo.95 = Estimate + qnorm(0.025)*Std..Error,
												ci.hi.95 = Estimate + qnorm(0.975)*Std..Error,
												ci.lo.90 = Estimate + qnorm(0.05)*Std..Error,
												ci.hi.90 = Estimate + qnorm(0.95)*Std..Error
)

treatment_tab$cov <- factor(c(0,1,0,1,0,1))
treatment_tab$Partisanship <- c("","","Republicans","Republicans","Democrats","Democrats")
treatment_tab$Partisanship <- factor(treatment_tab$Partisanship,levels=c("","Republicans","Democrats"),ordered=T)
treatment_tab$Subgroups <- c("Full Sample","Full Sample","Partisan\nSubgroups","Partisan\nSubgroups","Partisan\nSubgroups","Partisan\nSubgroups")

treatment_tab

## Figure 5:
plot_treat_nimby_intx <- ggplot(treatment_tab, aes(x=Partisanship,group=cov)) +
	geom_errorbar(aes(ymin=ci.lo.95, ymax=ci.hi.95),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=0.35) +
	geom_errorbar(aes(ymin=ci.lo.90, ymax=ci.hi.90),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=1.5) +
	geom_point(aes(y=Estimate,shape=cov),position=position_dodge(width=0.2),
						 data=treatment_tab,size=4) + 
	geom_hline(aes(yintercept=0),lty=2) + 
	# geom_vline(aes(xintercept=1.5),lty=1,lwd=3) + 
	facet_grid("Subgroups",scales = "free_x",space = "free_x") + 
	# xlab(NULL) +
	ylab("Effect of spatial self-interest\non policy support") +
	theme_bw() + 
	# scale_y_continuous(limits=c(-0.175, 0.005)) + 
	scale_x_discrete("") +
	scale_shape_discrete("Covariate-\nadjusted",labels=c("No","Yes")) + 
	# coord_flip() + 
	theme(text=element_text(colour=black,
													# family="CM Roman",
													size=15),
				legend.position = "right")

## Output:
fname = "plot_treat_nimby_intx.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_treat_nimby_intx, width=9, height=6)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=9, height=6)
	print(plot_treat_nimby_intx)
	dev.off()
}


##### Treatment effects by other moderators

treatment_tab <- data.frame(
	rbind(
		c(summary(nimby_cons_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_cons_nocov))),
		c(summary(nimby_cons)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_cons))),
		c(summary(nimby_lib_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_lib_nocov))),
		c(summary(nimby_lib)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_lib))),
		
		c(summary(nimby_white_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_white_nocov))),
		c(summary(nimby_white)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_white))),
		c(summary(nimby_black_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_black_nocov))),
		c(summary(nimby_black)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_black))),
		
		c(summary(nimby_hiinc_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_hiinc_nocov))),
		c(summary(nimby_hiinc)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_hiinc))),
		c(summary(nimby_loinc_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_loinc_nocov))),
		c(summary(nimby_loinc)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_loinc))),
		
		c(summary(nimby_hineed_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_hineed_nocov))),
		c(summary(nimby_hineed)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_hineed))),
		c(summary(nimby_loneed_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_loneed_nocov))),
		c(summary(nimby_loneed)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_loneed))),
		
		c(summary(nimby_pers1_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_pers1_nocov))),
		c(summary(nimby_pers1)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_pers1))),
		c(summary(nimby_pers0_nocov)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_pers0_nocov))),
		c(summary(nimby_pers0)$coefficients["distanceDistance:\nNear",c("Estimate","Std. Error")],n=nrow(model.frame(nimby_pers0)))
	)
)
treatment_tab <- mutate(treatment_tab,
												ci.lo.95 = Estimate + qnorm(0.025)*Std..Error,
												ci.hi.95 = Estimate + qnorm(0.975)*Std..Error,
												ci.lo.90 = Estimate + qnorm(0.05)*Std..Error,
												ci.hi.90 = Estimate + qnorm(0.95)*Std..Error
)

treatment_tab$cov <- factor(c(0,1,0,1,
															0,1,0,1,
															0,1,0,1,
															0,1,0,1,
															0,1,0,1))
treatment_tab$level <- c("Conservatives","Conservatives","Liberals","Liberals",
												 "White","White","Black","Black",
												 "Above-Median","Above-Median","Below-Median","Below-Median",
												 "Above-Median","Above-Median","Below-Median","Below-Median",
												 "Yes","Yes","No","No"
)
treatment_tab$level <- factor(treatment_tab$level,levels=c("Conservatives","Liberals","White","Black","Above-Median","Below-Median","Yes","No"),ordered=T)
treatment_tab$Subgroups <- c("Ideological\nSubgroups","Ideological\nSubgroups","Ideological\nSubgroups","Ideological\nSubgroups",
														 "Racial\nSubgroups","Racial\nSubgroups","Racial\nSubgroups","Racial\nSubgroups",
														 "Income\nSubgroups","Income\nSubgroups","Income\nSubgroups","Income\nSubgroups",
														 "Overdose-Rate\nSubgroups","Overdose-Rate\nSubgroups","Overdose-Rate\nSubgroups","Overdose-Rate\nSubgroups",
														 "Know Someone\nw/ Addiction","Know Someone\nw/ Addiction","Know Someone\nw/ Addiction","Know Someone\nw/ Addiction")

treatment_tab

## Figure A3:
plot_treat_nimby_intx_app <- ggplot(treatment_tab, aes(x=level,group=cov)) +
	geom_errorbar(aes(ymin=ci.lo.95, ymax=ci.hi.95),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=0.35) +
	geom_errorbar(aes(ymin=ci.lo.90, ymax=ci.hi.90),position=position_dodge(width=0.2),
								data=treatment_tab,
								width=.0, size=1.5) +
	geom_point(aes(y=Estimate,shape=cov),position=position_dodge(width=0.2),
						 data=treatment_tab,size=4) + 
	geom_hline(aes(yintercept=0),lty=2) + 
	# geom_vline(aes(xintercept=1.5),lty=1,lwd=3) + 
	facet_grid("Subgroups",scales = "free_x",space = "free_x") + 
	# xlab(NULL) +
	ylab("Effect of spatial self-interest\non policy support") +
	theme_bw() + 
	# scale_y_continuous(limits=c(-0.175, 0.005)) + 
	scale_x_discrete("") +
	scale_shape_discrete("Covariate-\nadjusted",labels=c("No","Yes")) + 
	# coord_flip() + 
	theme(text=element_text(colour=black,
													# family="CM Roman",
													size=15),
				legend.position = "right")

## Output:
fname = "plot_treat_nimby_intx_app.pdf"
if (Sys.info()['sysname'] == 'Windows'){
	ggsave(fname, plot_treat_nimby_intx_app, width=15, height=6)
	embed_fonts(fname, outfile=fname)
} else {
	CairoPDF(fname, width=15, height=6)
	print(plot_treat_nimby_intx_app)
	dev.off()
}


## Subgroup means:
dist_treat_all <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear"],y=a$dist_support_bi[a$distance=="Distance:\nFar"])

dist_treat_hiinc <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_income==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_income==1])
dist_treat_loinc <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_income==0],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_income==0])
dist_treat_hideath <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_death_rate==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_death_rate==1])
dist_treat_lodeath <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_death_rate==0],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_death_rate==0])
dist_treat_dem <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$dem==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$dem==1])
dist_treat_rep <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$rep==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$rep==1])

dist_treat_pers1 <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$personal_bi==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$personal_bi==1])
dist_treat_pers0 <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$personal_bi==0],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$personal_bi==0])

dist_treat_hiinc_dem <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_income==1 & a$dem==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_income==1 & a$dem==1])
dist_treat_loinc_dem <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_income==0 & a$dem==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_income==0 & a$dem==1])
dist_treat_hiinc_rep <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_income==1 & a$rep==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_income==1 & a$rep==1])
dist_treat_loinc_rep <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_income==0 & a$rep==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_income==0 & a$rep==1])

dist_treat_hideath_dem <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_death_rate==1 & a$dem==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_death_rate==1 & a$dem==1])
dist_treat_lodeath_dem <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_death_rate==0 & a$dem==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_death_rate==0 & a$dem==1])
dist_treat_hideath_rep <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_death_rate==1 & a$rep==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_death_rate==1 & a$rep==1])
dist_treat_lodeath_rep <- t.test(x=a$dist_support_bi[a$distance=="Distance:\nNear" & a$above_death_rate==0 & a$rep==1],y=a$dist_support_bi[a$distance=="Distance:\nFar" & a$above_death_rate==0 & a$rep==1])


table_dist_diffs <- rbind(
	cbind("All Respondents",round(dist_treat_all$estimate[2],3),round(dist_treat_all$estimate[1],3),round(dist_treat_all$estimate[1]-dist_treat_all$estimate[2],3),round(dist_treat_all$p.value,3)),
	cbind("","","",paste0("(",round(dist_treat_all$conf.int[1],3),", ",round(dist_treat_all$conf.int[2],3),")"),""),
	
	cbind("n",nrow(a[which(a$distance=="Distance:\nFar"),]),nrow(a[which(a$distance=="Distance:\nNear"),]),"",""),
	cbind("","","","",""),
	cbind("Above Median Income",round(dist_treat_hiinc$estimate[2],3),round(dist_treat_hiinc$estimate[1],3),round(dist_treat_hiinc$estimate[1]-dist_treat_hiinc$estimate[2],3),round(dist_treat_hiinc$p.value,3)),
	cbind("","","",paste0("(",round(dist_treat_hiinc$conf.int[1],3),", ",round(dist_treat_hiinc$conf.int[2],3),")"),""),
	
	cbind("n",nrow(a[which(a$distance=="Distance:\nFar" & a$above_income==1),]),nrow(a[which(a$distance=="Distance:\nNear" & a$above_income==1),]),"",""),
	cbind("","","","",""),
	
	cbind("Below Median Income",round(dist_treat_loinc$estimate[2],3),round(dist_treat_loinc$estimate[1],3),round(dist_treat_loinc$estimate[1]-dist_treat_loinc$estimate[2],3),round(dist_treat_loinc$p.value,3)),
	cbind("","","",paste0("(",round(dist_treat_loinc$conf.int[1],3),", ",round(dist_treat_loinc$conf.int[2],3),")"),""),
	
	cbind("n",nrow(a[which(a$distance=="Distance:\nFar" & a$above_income==0),]),nrow(a[which(a$distance=="Distance:\nNear" & a$above_income==0),]),"",""),
	cbind("","","","",""),
	
	
	cbind("Above Median Overdose Rate",round(dist_treat_hideath$estimate[2],3),round(dist_treat_hideath$estimate[1],3),round(dist_treat_hideath$estimate[1]-dist_treat_hideath$estimate[2],3),round(dist_treat_hideath$p.value,3)),
	cbind("","","",paste0("(",round(dist_treat_hideath$conf.int[1],3),", ",round(dist_treat_hideath$conf.int[2],3),")"),""),
	
	cbind("n",nrow(a[which(a$distance=="Distance:\nFar" & a$above_death_rate==1),]),nrow(a[which(a$distance=="Distance:\nNear" & a$above_death_rate==1),]),"",""),
	cbind("","","","",""),
	cbind("Below Median Overdose Rate",round(dist_treat_lodeath$estimate[2],3),round(dist_treat_lodeath$estimate[1],3),round(dist_treat_lodeath$estimate[1]-dist_treat_lodeath$estimate[2],3),round(dist_treat_lodeath$p.value,3)),
	cbind("","","",paste0("(",round(dist_treat_lodeath$conf.int[1],3),", ",round(dist_treat_lodeath$conf.int[2],3),")"),""),
	
	cbind("n",nrow(a[which(a$distance=="Distance:\nFar" & a$above_death_rate==0),]),nrow(a[which(a$distance=="Distance:\nNear" & a$above_death_rate==0),]),"",""),
	cbind("","","","",""),
	
	cbind("Democratic Respondents",round(dist_treat_dem$estimate[2],3),round(dist_treat_dem$estimate[1],3),round(dist_treat_dem$estimate[1]-dist_treat_dem$estimate[2],3),round(dist_treat_dem$p.value,3)),
	cbind("","","",paste0("(",round(dist_treat_dem$conf.int[1],3),", ",round(dist_treat_dem$conf.int[2],3),")"),""),
	
	cbind("n",nrow(a[which(a$distance=="Distance:\nFar" & a$dem==1),]),nrow(a[which(a$distance=="Distance:\nNear" & a$dem==1),]),"",""),
	cbind("","","","",""),
	cbind("Republican Respondents",round(dist_treat_rep$estimate[2],3),round(dist_treat_rep$estimate[1],3),round(dist_treat_rep$estimate[1]-dist_treat_rep$estimate[2],3),round(dist_treat_rep$p.value,3)),
	cbind("","","",paste0("(",round(dist_treat_rep$conf.int[1],3),", ",round(dist_treat_rep$conf.int[2],3),")"),""),
	
	cbind("n",nrow(a[which(a$distance=="Distance:\nFar" & a$rep==1),]),nrow(a[which(a$distance=="Distance:\nNear" & a$rep==1),]),"",""),
	cbind("","","","",""),
	
	cbind("Know Someone with Addiction",round(dist_treat_pers1$estimate[2],3),round(dist_treat_pers1$estimate[1],3),round(dist_treat_pers1$estimate[1]-dist_treat_pers1$estimate[2],3),round(dist_treat_pers1$p.value,3)),
	cbind("","","",paste0("(",round(dist_treat_pers1$conf.int[1],3),", ",round(dist_treat_pers1$conf.int[2],3),")"),""),
	
	cbind("n",nrow(a[which(a$distance=="Distance:\nFar" & a$personal_bi==1),]),nrow(a[which(a$distance=="Distance:\nNear" & a$personal_bi==1),]),"",""),
	cbind("","","","",""),
	cbind("Don't Know Someone with Addiction",round(dist_treat_pers0$estimate[2],3),round(dist_treat_pers0$estimate[1],3),round(dist_treat_pers0$estimate[1]-dist_treat_pers0$estimate[2],3),round(dist_treat_pers0$p.value,3)),
	cbind("","","",paste0("(",round(dist_treat_pers0$conf.int[1],3),", ",round(dist_treat_pers0$conf.int[2],3),")"),""),
	
	cbind("n",nrow(a[which(a$distance=="Distance:\nFar" & a$personal_bi==0),]),nrow(a[which(a$distance=="Distance:\nNear" & a$personal_bi==0),]),"",""),
	cbind("","","","","")
	
)
colnames(table_dist_diffs) <- c("Subset","Mean[Far]","Mean[Near]","Treatment effect (CI)","p-value of difference")

## Output Table A8:
print(xtable(table_dist_diffs),floating=F,file = "table_dist_diffs.tex",include.rownames=F,booktabs=T)

